%**************************************************************************
% DISE model
% Aneli Bongers and Jose L. Torres (2025)
% Paper: "Optimal path for orbital debris"
% File: dynamics.m
% Function dynamics
%**************************************************************************
function [f,varargout]  =   dynamics(x,v,i,Params,varargin);
nouextra=max(nargout)-1;
ninextra=max(nargin)-4;
%**************************************************************************
% States
%**************************************************************************
k       =   x(1);
s       =   x(2);
W       =   x(3);
Z       =   x(4);
F       =   x(5);
Welfare =   x(6);
%**************************************************************************
% Endogenous solutions
%*************************************************************************
srk=v(1);
srs=v(2);
%**************************************************************************
%Calibration
%**************************************************************************
%Params=Space1999;
%**************************************************************************
%Variables
%**************************************************************************
y=Params.a(i)*k^Params.alpha1*s^Params.alpha2*(Params.N(i))^(1-Params.alpha1-Params.alpha2);
cc=(1-srk-srs)*y;
S=Params.mu*s/Params.q(i);
D=W+Z+F;
destroyed=Params.theta*D*s;
l=Params.m(i)*srs*y;
L=Params.mu*(1-Params.m(i))*srs*y/Params.eta;
XX=Params.theta*D*S;

 switch ninextra
     case 1
         u=varargin{1};
         c=u(1);
     case 0
         c=cc;
 end

k_next=(1-Params.deltak)*k+srk*y;
s_next=(1-Params.deltas)*s+Params.q(i)*(1-Params.m(i))*srs*y-destroyed;
W_next=(1-Params.deltaw-Params.epsilonw)*W-Params.theta*(D+S)*W+Params.xi*Params.deltas*S;
Z_next=(1-Params.deltaz-Params.epsilonz)*Z-Params.theta*(D+S)*Z+Params.psi*L;
F_next=(1-Params.deltaf)*F+(1+Params.nu1)*Params.phiw*Params.epsilonw*W+(1+Params.nu1)*Params.phiz*Params.epsilonz*Z+(1+Params.nu1)*Params.omega*L+(1+Params.nu1)*Params.gammas*XX+(1+Params.nu1)*Params.gammaw*Params.theta*D*W+(1+Params.nu1)*Params.gammaz*Params.theta*D*Z;
Welfare_next=Welfare-Params.N(i)*((((c/Params.N(i))^(1-Params.sigma)-1)/(1-Params.sigma)))/(1+Params.rho)^(i-1);

f=[k_next;s_next;W_next;Z_next;F_next;Welfare_next];

if nouextra == 1
     varargout{1}=cc;
end
